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METHO DS OF DIGITAL FILTERING ANU MULTI-DIME MSIO NAL DATA COMPRESSION USING THE FAREY 
t?07a3KXTTOE'7QsJD ARITHMETIC, FAN, AND MODULARWAVhEItS ~" ~~" 



It is well-known that digital filtering of digitally represented samples of analog data can 
be simplified and improved by performing the required processing in some transformed 
domain. Using^both a transform and its inverse to simplify the digital or other manip- 
ulation of data, one first transforms the data from its initial form to an intermediate 
form, then applies a specified manipulation, and finally applies the inverse transform to 
recover the data in its initial form from its intermediate form. By way of example and 
not by way of limitation: the initial form of the data may be digital samples of ana- 
log data, and the intermediate form of the data may be the coefficients of the Fourier 
transform; the initial form of the data may he the coefficients of the Fourier transform, 
and the intermediate form of the data may be the coefficients in some wavelet or other 
expansion; the initial form of the data may be an analog signal, and the intermediate 
form of the data may be some digital representation of it; the initial form of the data 
may be digital samples of an analog signal, and the intermediate form of the data may 
be a quantization of the digital samples. For instEince, transform coding methods for 
data compression in signal or data processing are realized in practice as transform of 
data/quantization, storage or transmission, de-quantization/reconstruction. Another 
class of examples consists of finite-impulse response digital filters where the input 
data is filtered using the indirect calculation of convolution given by Fourier trans- 
form/multiplication/inverse Fourier transform. For analog data, sampfing and filtering 
together determine fidelity and speed of manipulation. 

The invention disclosed here provides a new method for sampling analog or digital 
input data which may be used to calculate new wavelet transforms as well as their 
inverse reconstruction algorithms. This immediately provides novel and efficient meth- . 
ods for data compression based on this new method of non-linear transform coding. 
In combination with these wavelet filters, the invention also provides new methods for 
calculating various classical transforms including the Fourier transform and its inverse. 
This immediately provides novel and efficient methods for digital filtering. By way 
of ex£imple and not of limitation, practical applications of the methods include non- 
speech audio compression, speech compression, speech recognition, speech synthesis, 
voice printing, audio filtering for hearing aids, still two-dimensional image compres- 
sion, moving video compression, video compression for purposes of telephony, precision 
Fourier analysis, precision trigonometry, denoising, interpolation, medical imaging, ge- 
ological imaging for recovery of oil or other resources, and other emission/detection 
apparatus. In specific applications, it may be necessary to store streams of input data 
in a buffer for separate processing or manipulation while concurrently collecting further 
input data whether analog or digital. It may also be useful in practice in particular 
applications to calculate and archive or otherwise store specific coefficients or constants 
and recover them from memory at run time rather than computing them at run time. 



Field of the Invention 



Description of Prior Art 



In the landmark paper "Orthonormal bases of compactly supported wavelets. Com- 
munications in Pure and Applied Mathematics 41 , pp. 909-996 (1988), Daubechies 
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presented a class of w :\ets which have been effectively applie' i numerous con- 
texts; see, for instance. US Patents 5,453,945 by Tucker et al., 5,605,575 by Williams, 
as well as Daubechies' monograph "Ten Lectures on Wavelets" SIAM (1992) and its 
references. In another landmark paper, "An algorithm for the machine calculation of 
complex Fourier series", Mathematics of Computation 19, pp. 297-301 (1965), J. W. 
Cooley and J. W. Tukey described a simplified algorithm for the calculation of Fourier 
series on the computer. The utility of this fast Fourier transform FFT is well known 
with thousands of important applications and implementations; see, for instance, US 
Patents 3,871,577 by Avellar et al., 3,881,100 by Works et al., 4,051,357 by Bonnerot, 
and "The DFT, an owner's manual for the discrete Fourier transform", SIAM (1995) 
by W. Briggs and Van Emden Henson and its references. A requirement for both the 
FFT method and for the multiscale filtering of the Daubechies' wavelets is that the 
sampling of input data must be equally spaced though further methods of adaptive 
refinement of grids have also been developed in each context. 

Other examples of transforms whose efficient digital implementations are important 
and to which the invention is relevant include the transforms of Hilbert, Haar, Laplace, 
Bessel, Laguerre, Hermite, Chebyshev, Hotelling, Mersenne, and Fermat; see, for in- 
stance, US Patents 3,891,443 by Lynch et al. and 4,093,994 by Nussbaumer. 



. ' ^ Invention Summary 

-■ f 

Three new families of wavelets, respectively called arithmetic wavelets, fan wavelets, and 
modular wavelets, are defined on the standard circle of radius one in the complex plane. 
Any complex-valued function defined on this circle may be written essentially uniquely 
as a complex linear combination of each family, and the approximation of specified input 
data as a finite linear combination for each family leads to a corresponding wavelet filter. 
Each filter depends upon a new quadrature on the circle, called the Farey quadrature, 
which relies on a novel method, of non-equally-spaced sampling. The output of each 
wavelet filter plays the role in the invention of an intermediate representation between 
input data and any one of a number of useful output transformations of data which 
may be computed from the intermediate quantity. 

With either analog or digital input data, it is convenient to think of the given in- 
put data as spatial data. The inverse wavelet transform or reconstruction algorithm 
requires calculating spatial data from wavelet data and admits an especially conve- 
nient implementation: the values of the spatial function at certain non-equally-spaced 
grids of arbitrarily fine spacing may be computed exactly using a specified finite set 
of wavelet coefficients. This reconstruction algorithm combines with the wavelet fil- 
ter into a binary cascade giving an efficient procedure for data compression which is 
especially well-suited to discontinuous input data and to iterative refinement of out- 
put data. Furthermore, the methods extend directly to procedures for compression 
of multi-dimensional data. Computer coding for the reconstruction algorithm itself 
is sufficiently abbreviated and low-level that it might be transmitted together with 
compressed data, for instance, for down-loading from satellite to home computer. 

For the calculation of classical transforms of spatial input data, it is convenient to think 
of the desired output data, for instance, the Fourier coefficients, as frequency data. The 
transform from spatial to frequency data is accoraphshed in two main steps, the first 
of which approximates the input data by wavelets using the wavelet filter. The second 
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step is the calculation 7ourier coefficients or other frequency dat . nerraSTJf fvaMetS 
and requires formulas derived in previous theoretical -work. Similarly, the calculation of 
the inverse transform depends upon the reconstruction algorithm and known formulas 
which express frequency data in terms of the wavelets. The calculation of a transform 
or its inverse admits an elegant implementation as a binary cascade. 



Brief Description of Figures 

FIG la illustrates the standard circle C of radius one in the complex plane, where 
C bounds the standard disk D of radius one. The complex numbers +1 and -1 
are also indicated, and the chord of C with these endpoints is labeled by the matrix 



FIG lb illustrates a chord in the top semi-circle of C labeled by a matrix ^4 = ^ " ^ ^ , 

together with two further "descendant" chords of C labeled by the matrix products 

UA and TA as indicated, %vhere J^=(j °^andT=Q J^. The three chords 

determine a triangle inscribed in C, and the vertices of this triangle are indicated as 
complex numbers. 

FIG Ic illustrates a chord in the bottom semi-circle of C labeled by a matrix A = 
(c d)' *°Sether with two further "descendant" chords of C labeled by the ma- 
trix products U~^A and T~'^A as indicated, where U~'^ = ^ l} ~ 

* "^^^ ^^'^^^ chords determine a triangle inscribed in C, and the vertices of 
this triangle are indicated as complex numbers. 

FIG 2 depicts a classical figure in mathematics called "the Farey tesselation in Klein's 
model of hyperboHc geometry". It is constructed beginning from the chord in Figure la 
by recursively taking descendants as in Figure lb and Figure Ic, respectively, on the 
top and bottom semi-circle of C. The figure itself is classical, but the sampling of 
input data on C at the endpoints of the chords of the Farey tesselation in increasing 
generation is a novel aspect of the invention disclosed here. 

FIG 3a illustrates the circle C, the chord labeled / and the two triangles in the Farey 
tesselation on either side of this chord. The vertices of these triangles are indicated 
inside the circle as complex numbers. Outside the circle are given four explicit com- 
binations of cosine and sine functions, one such function next to each circular arc 
determined by the vertices, where & is the usual angular coordinate on the circle. A 
function "QiiGi) is uniquely determined by the figure, where on each such circular arc, 
"di^B) takes the values of the neeirby combination of cosine and sine functions. The 
other parts Figures 3b-3e likewise determine analogous functions "Qa (^) for other ma- 
trices A = These functions taken together form the system of wavelets on 
which the invention is based, and they are unique to this invention. 
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FIG 3b illustrates th- .lord in the top semi-circle of C labeled th the matrix A, 
where a > c. As before, the figure determines a function on the circle denoted t9,4(5). 

FIG 3c illustrates the chord in the top semi-circie of C labeled with the matrix A, 
where c >= a. As before, the figure determines a function on the circle denoted 

FIG 3d illustrates the chord in the bottom semi-circle of C labeled with the matrix A, 
where a > —c. As before, the figure determines a function on the circle denoted -OaIQ). 

FIG 3e illustrates the chord in the bottom semi-circle of C labeled with the matrix 
A where — c >= a. As before, the figure determines a function on the circle denoted 

FIG 4 illustrates the notation near the chord labeled by the matrix A = ^) 
the top semi-circle of C; this notation will be useful in several subsequent discussions. 

FIG 5 provides a flow chart for the main driving routine in calculating the Fourier 
transform. 

FIG 6 provides a flow chart for the main recursive routine in calculating the Fourier 
transform. 

PIG 7 provides a flow chart for the procedure of updating coefficients in calculating 
the Fourier transform. 

FIG 8 provides a flow chart for the procedure of processing terminal arrows of the 
recursion in calculating the Fourier transform. 

FIG 9 provides a flow chart for the main driving routine in calculating the inverse 
Fourier transform. 

FIG 10 provides a flow chart for the main recursive routine in calculating the inverse 
Fourier transform. 

FIG 11 illustrates the region U consisting of all complex numbers with positive imag- 
inary part. Map the disk D toU via the function z ^ i{z + l)/(z — 1); if two points in 
C are the endpoints of a chord in Klein's model of the Farej' tesselation as in Figure 2, 
then draw the semi-circle in U passing through the corresponding two points which is 
perpendicular to the real axis. This produces another classical figure in mathematics, 
called the "Farey tesselation in the upper half-space model of hyperbolic geometry", 
which is indicated in Figure 11. 



Description of Preferred Embodiments 

In order to disclose the methods, it is necessary to develop some notation and termi- 
nology. To this end, consider the standard disk D of radius one about the origin in 
the complex plane with its usual complex coordinate z = x + iy, where i — and 
x,y 3xe real numbers, together with its boundary circle C as illustrated in Figure la. 
An arrovj is by definition an oriented chord of the circle C labeled by a two-by-two 

matrix f ^ ^ j , where the initial point of the chord is given by the complex number 
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f^g e C, and the final >int of the chord is given by 6 C. A 2cial arrow called 
the doe (short for "distinguished oriented edge") is given by the diameter of C with 

endpoints — labeled by the matrix = ^ Figure la. Also required 

subsequently are the matrices S = (j^ ^ and [/" = ^ ^ ° j , = J 1 ) ' 

any integer n, which satisfy the identities 5^ = -/, SUS = -T"^, STS = -U''^ of 
matrix products. 

Inductively define two separate binary cascades, one of top arithmetic arrows and an- 
other of bottom arithmetic arrows, as follows: 

Basis: The doe, which is illustrated in Figure la, is both a top arithmetic 
arrow and a bottom arithmetic arrow. 

Top Induction: As in Figure lb, given the top arithmetic arrow labeled by 
~ c d ) ' ^^^^ arithmetic arrows labeled by 

- UA^( I ,_^^') and T^=f^ + ^ ■ 
\a + c b + dj \ c d J 

Bottom Induction: As in Figure Ic, given the bottom arithmetic arrow 
labeled by ^ = ^ ^ ^ ^ , there are also bottom arithmetic arrows labeled 
by 



U-'A=( ° and r-^^=f°- 

\c — ad—bj \ c 



It is a classical theorem in mathematics that the chords underlying any two arithmetic 
arrows meet only at their endpoints, if at all, and that the family of all such chords 
decomposes the disk D into an infinite collection of triangles with vertices in the circle 
C as is illustrated in Figure 2; this figure is called "the Farey tesselation in Klein's 
model of hyperbolic geometry" . There is no suitable reference for this procedure in the 
literature, and it is for completeness that the explicit recursive definition of arithmetic 
arrows has been given here. An arithmetic arrow is said to be of generation g if it arises 
from g apphcations of the inductive clause starting from the doe in case of either top 
or bottom arrows; the doe heis generation zero by convention. 

Consider the infinite set Q C C consisting of the endpoints of chords underlying all 
arithmetic arrows, and say a point of Q is of generation g if it is the endpoint of the 
chord underlying an arrow of generation g and no less; —1 and +1 are of generation 
zero by convention. There is a canonical enumeration of the generation g points of Q 
clockwise in C starting from —1 and the associated enumeration of Q itself by increasing 
generation. For example, the ordering begins with: 

-1, +1, +i, 

Given input data on the circle C, the sampling of it at the points Q C. C in this order 
of increasing generation will be referred to as the Farey quadrature. In practice, one 
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might interpolate given ta as required at the sampling points of tl "^arey quadrature 
or restrict the Farey quadrature to a circular arc in C. (In the subsequent discussion 
of mathematical basis, the sample points of the Farey quadrature will be seen to corre- 
spond to the rational points in the real line enumerated in increasing order, where the 
generation is the length of a corresponding continued fraction expansion.) 

Notice that the doe separates the triangle with vertices — 1,4-1, -i from the triangle 
with vertices —1, +1, +i. In the same way, an arbitrary arithmetic arrow separates two 
triangles complementary to all the arithmetic arrows in D. The endpoints of chords 
of these triangles in various cases are illustrated in Figures 3a-3e as functions of the 

matrix A = labeling the arithmetic arrow. Figure 3a depicts the doe itself 

with A — Iy Figure 3b and 3c depict top arrows with a> c and c > a, respectively, and 
Figure 3d and 3e depict bottom arrows with a > -c and -c > a,' respectively. The 
notation inside the circles indicates the endpoints of the chords as complex numbers, 
and the notation outside the circles in each case of Figures 3a-3e is yet to be explained. 

To this end, define a trigonometric function to be a real-valued function 

f{6) = acose+ Psine + j, 

where 6 is the usual angular coordinate on the unit circle C and a, /5, 7 are real numbers. 
Next, define a piecewise trigonometric function on the circle to be a function f{6) so 
that there is a finite collection of circular arcs Ij C C, for j = 1, . . . , J, any two of which 
meet at common endpoints in C, if at all, together with a collection of trigonometric 
functions tj{e) for i = 1,. .., J, so that f{d) = ij{6) if 6 € Ij. In other words, C is 
decomposed into a finite number of non-overlapping circular arcs, and on each such 
arc, / takes the values of a trigonometric function. 

There is a particular family of piecewise trigonometric functions, one such arithmetic 
wavelet ^/i(^) for each matrix A labeling an arithmetic arrow. These functions are 
defined in Figures 3a-3e, where the endpoints of the circular arcs are as described 
before, and each trigonometric function is written outside the circle next to its corre- 
sponding circular arc. There are thus five cases in the definition of arithmetic wavelet 
corresponding to the five cases in Figures 3a-3e already discussed, and in each case, the 
arithmetic wavelet itself is a piecewise trigonometric function defined using exactly four 
circular arcs. Arithmetic wavelets are a fundamental new ingredient of this invention 
and no doubt seem entirely ad hoc axid unmotivated at this juncture of the discussion. 
The mathematical basis for them will be given later, but it is worth pointing out now 
that each arithmetic wavelet '&a(0) is a once-continuously differentiable function taking 
value zero at the points —1, +1, —i of the circle C, that is, for 6 = 0, — 7r/2, — tt, except 
for ■du-i{6) and which take the value zero only at the points ±1. 

In order to explicate the definition of arithmetic wavelets, notice that any two-by-two 



matrix A 




determines a self-mapping Ma : C — > C of the circle C as follows: 



_ {ta -F fe) - i{tc -h d) 
^\t + i) ~ ita + b) + i(tc + dy 

and in fact, the image of the endpoints, —1 and -hi, of the doe under Ma are the 
endpoints of the arithmetic arrow with label A. Construct from the function di{9) 

6 
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defined on C in Figure I .he corresponding vector field on C given di{9) ■§g, where 
^ is the usual unit tangent vector field to C. 

There is an explicit formula to be given presently for the transformation of suitable 
vector fields under the change of coordinates on C given by Ma- Namely, the vector 
field t{e)-§Q = [a cos e+jSsm ^'+7]^ transforms under Ma to the vector field t^{6)-§g, 
where 

^ [2(cd - ab)0 + {d^ - b'^Xa + 7) + (a^ - c^) - 7)] cos 6 
(e) = I + [{ad + bc)p + bd{a + 7) - ac{a - 7)] sin ^ | . 

+ i {2{cd + ab)0 + (d^ + b^){a + 7) - (o^ + c^){a - 7)] 

More generally, given a piecewise trigonometric function j{9) taking the values of the 
trigonometric function 10) on the circular arc Ij as before, the corresponding vector 
field j{Q)^ transforms under Ma to the vector field which on the circular arc Ma(/j) 
isgivenbyt/{0)^. 

Changing coordinates on C by Ma using this formula transforms the vector field 
^/(^) ^ to another vector field ■Qa{&)^, thereby defining the function on C. 

Finally except for A = and A = the corresponding arithmetic wavelet is 
given by 

^a(^) = <iA{9) - [a cos ^ + /3 sin Q + 7], 

where o;,;3,7 are chosen so that i^(7r) = i^CO) = d{-Tt/2) = 0; the two special cases are 

du-i(0)='du-i{e) + 2sme and = (e) - 2 sin ^ 

(which are defined in this manner to guarantee a symmetry of the upper and lower 
semi-circles of C). This is the abstract definition of arithmetic wavelets which is made 
explicit in Figure 3. 

It is furthermore convenient to define 

(dAiey, if Ay^U-\T-^, 

[^T-^{9); ifA = T-\ 
so that for every label A on an arithmetic arrow, t?A{7r) = ^/i(0) = -i5^(-7r/2) = 0. 
The arithmetic wavelets enjoy an important renormalization property, namely, 
^BA{MAie)) = MM^Bie), 

where the prime denotes the usual derivative and where it is required that either both 
matrices B and A are matrix products of {7"*"Sr+^ or both matrices B and A are 
matrix products of U~^,T~^. 

Just as a function / = f{6) may be expressed in its representation as a Fourier series / ~ 
Yin Cn^'"^, so too it may be expressed as a series f ~J2 ^a'^a (9) of arithmetic wavelets, 
which together form a linearly independent set. The sum is over all arithmetic arrows, 
and the coefficients ca are called the arithmetic wavelet coefficients. The calculation 
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of arithmetic wavelet c acients from the input data / is called th rithmetic wavelet 
filter. The arithmetic wavelet coefficients are not quite uniquely determined by / as 
will be explained. 

The reconstruction algorithm or inverse wavelet transform provides a method for cal- 
culating function values at points of the circle from arithmetic wavelet coefficients. To 
determine the value of the sum X) ^a'Oa at a specified point g e Q Q C, draw a line 
segment from g to the origin in D; this segment is contained in a finite collection of tri- 
angles complementary to all arithmetic arrows, whose boundary edges taken together 
form a finite set of chords; the corresponding finite set of arithmetic arrows A enu- 
merates the only coefficients which affect the value of X^^^^ sa^a at g. This is an 
important finiteness property of arithmetic wavelets which is crucial to the efficiency of 
the reconstruction algorithm. 

This completes the general discussion of arithmetic wavelets, the first family of new 
wavelets disclosed here, and the two further families of new wavelets are next defined. 

For each label A on an arithmetic arrow, there are two fan wavelets 

n>0 
n>0 

as well as two modular wavelets 

n>0 

i^sA{e) = n Ir-^Aie). 

n>0 

It follows immediately from the definitions that there are the following relations among 
the wavelets for any label A on an arithmetic arrow. 

<t>A{6) = ^u-^A{e)-MO), 4>sa{0) = -^u-^sAm-i^sAie), 
^A{d) = Me) - <i>uAie) = 4>sA{e)-4>usA{6) 

= ■>Pu-^a{0)-2M6)+iPua{6) = '4>u-^sA{&)-'2i^sA{d)+i'usA{0)- 

On the other hand as is not at all obvious from the definitions, 4>i and ipi are given by 
the following explicit formulas as piecewise trigonometric functions 

{2 sin ^; if 0 < ^ < tt, 

-2 cos 0 - 2 sin 6' - 2; if tt < < 37r/2, 
-2 cos ^ -f 2 sin 9 + 2; if 37r/2 <9<2'k, 

and 

, , f 2 - 2 cos ^; if 0 < e < TT, 
^^(^)- = \0; if7r<^<27r. 

In fact, modular wavelets arise from il)! in exactly the same manner that arithmetic 
wavelets arise from •^/, namely, transforms under Ma to the vector field 
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'P*Ai^)w-' "^Aid) = i(^) — [oc COS 6 + 13 sin 5 + 7], where a, 7 are chosen to 
guarantee that ^a(O) = VJa{-7^) = VA(-7i'/2) = 0: transforming ipi{6)-§Q using the 
change of coordinates given by Ms a likewise gives rise to tpsA(d) ■ The fan wavelets 
also arise from <pj in the analogous manner. 

In contrast to arithmetic wavelets, neither the fan nor modular wavelets are linearly 
independent, and in each case, there is one linear relation for each arithmetic arrow, 
namely, 

<pA - 4>UA = 4'SA - ^USAi 
■>PU-^A - ^■^A + i>UA = i>U-^SA - ^i^SA + tpUSA- 

On the other hand, the following collection of modular wavelets forms a basis 
{tpA : A labels a top arithmetic arrow and c > 0} 
U {tpSA ■ ^ labels a top arithemtic arrow and a > c} 
U {4?A '• A labels a bottom arithemtic arrow and — c > a} 
U {fpSA ' A labels a bottom arithemtic arrow and a > — c); 

there is an analogous basis for the fan wavelets. 

Both fan and modular wavelets enjoy precisely the same finiteness property as arith- 
metic wavelets. Furthermore, both fan and modular wavelets enjoy renormalization 
properties, namely, 

4>BAiMAm = MM Me), 
^ba{Ma{9)) = M'A{e)Me), 

where the requirements on Ay B are as before. 

This completes the necessary abstract definitions. The methods will first be presented 
using arithmetic wavelets for real-valued input functions with extensions to modular 
and fan wavelets and complex- valued or multi-dimensional input/output data to be 
described later. As a point of notation throughout these discussions, there will be the 
tacit assumption that the entries in the matrix denoted A are given by 



The Fourier Transform 

Included in the section entitled "Source Codes" is an implementation in the computer 
language C of a preferred embodiment of the method described in this section employing 
arithmetic wavelets. 

Suppose that f(9) is some specified real-valued input function defined on the circle. 
The calculation of the Fourier transform of / proceeds by two main steps: 

Stepl : Choose a finite set S of arithmetic arrows and approximate / w 
X^yies ^A^A calculating the coefficients ex in a manner to be described in 
terms of the values of / at specified points using the Farey quadrature. 

Step 2: Substitute into the expression in Step 1 the known Fourier expan- 
sions ^a{9) ~ Y^n e'"^ in order to derive the approximation , 

Aes 
9 
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The coefficients require n Step 2 are known theoretically and gi\ for n 0, ±1 by 



Furthermore, the coefficients c^,cf ,cfi are given by 

TTC^i = Ohic^ -d^± 2icd) + e^lib - df - (a - cf T 2i(6 - ^)(a - c)]/2 
+ ^t(a2 - ± 2m6) + ee[{b + df - {a + cf ^^2i{b + d){a + c)]/2, 

1^4 = 26h{c^ + d"^) - Griib-df + ia-cf] 
+ 20t(a2 + 62) _ 9e[{b + df + {a + cfl 

■where the angles are 



6, - tan-i (-^^tt^^^fl±SL.\ 



,, = tan-( .. = tan--( /^-/)("-) ). 

It may be useful in practice in particular applications to calculate and archive or oth- 
erwise Store the constants required in the previous expressions for and recover them 
from memory at run time rather than computing them at run time. 

The Fourier coefficients co,c±i of / thus require special treatment which does not 
present additional challenges. For simplicity, suppose that the input function f{9) 
has these three coefficients fixed by the condition that f{6) = 0 for ^ = 0, — 7r/2, — tt. 
More precisely, if S{9) is any real-valued function on the circle, then let f{9) denote its 
normalization defined by f {0) - f{6) + a cos ^ -F /? sin ^ -j- 7, where a, ^, 7 are chosen 
so that / takes value zero at 0, — 5r/2, — The specified simplifying assumption on / 
thus amounts to the condition that f — J. 

The procedure for calculating wavelet coefficients is slightly different depending upon 
whether the corresponding arithmetic arrow is a top or bottom arrow. However, the 
antipodal map 6 t-^ -k + 6 transforms the bottom of the circle C to the top in such a 
way that it suffices to discuss only the case of top arrows. More precisely, the algorithm 
for bottom arrows arises from that for top arrows by: 

• replacing the input data value f(6) with the input data value f(0 + tt), 

• replacing the matrix A in the described manipulations for top arrows 
by the conjugate matrix SAS for bottom arrows, and 

• replacing the trigonometric function acos 0 + 0sin d + j for top arrows 
by the trigonometric function —acos 6 - 0sin ^ -f- 7 for bottom arrows. 
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The algorithm and cor. ' iter coding for bottom arrows is thereft derived without 
difficulty from that for top arrows, and the subsequent discussion is specialized without 
loss of generality to top arrows. 

In addition to the function f{6) defined on the circle, suppose that there is also specified 
a bandwidth of interest N (i.e. the method should calculate the Fourier coefficients Cn 
only in the range jr?! < N) and a tolerance EPS (where contributions to the specified 
Fourier coefficients are regarded as negligible if they are of magnitude less than EPS). 
There are numerous variants of the method which depend upon other parameters which 
may be tailored to a given application, and some these variants will be discussed further. 

Steps 1 and 2 are merged into a single binary cascade which keeps track of an ongoing 
approximation to the Fourier coefficients as follows, A cascade element is called an 
arrow-structure and is defined to consist of the specification of 

• an arithmetic arrow A of generation 

• the wavelet coefficient ca on "Oa in the approximation of /, 

• the coefficients a, /?, 7 of a trigonometric function defined as follows. The end- 
points of the chord corresponding to A decompose C into two circular arcs, ex- 
actly one of which contains none of the points -1, +1, -i. The sum YIbt^a ^b^b 
over all arithmetic arrows B of generation at most g except for A itself is given 
on this interval by the particular trigonometric function a cos 0 + /? sin ^ 4- 7. 

The coefficients a, (3, 7 in an arrow-structure keep a lagged/updated running tally of 
the overall effect of what has come before it in the cascade; as a result of this technique, 
the method requires essentially no memory other than the storage of the running ap- 
proximation to Fourier coefficients and the stack required for the recursive computation 
in the cascade of arrow-structures. 

The next technical point involves a regularization in the calculation of arithmetic 
wavelet coefficients; this regularization is required because the approximation to / as a 
linear combination X^^^^ Ca^x is not uniquely determined. Given an arrow-structure 
in the specified notation and referring to Figure 4, when sampling the input data at 
the point in the Farey enumeration with corresponding angular coordi- 

nate do, the one real input datum f{do) must determine the two coefficients euA and 
cta- According to the formulas presented in Figures 3a-3e, these new coefficients are 
constrained by 

4. ^i^UA - erx) 8 Ca sgn (a - c) 

(6-Hd)2-F(a-hc)2 {b + d)-^ + {a + cY' 

where sgn (a — c) is a sign defined to be +1 if a > c and -1 otherwise. There remains, 
however, one real degree of freedom in the specification of euA, cta, and it is eliminated 
by demanding the best least-squares fit to the values of / at the next generation of points 
(2&+d)-i(2a+c) ' (bt2d)-i(a+2c) '^'^^^ respective angular coordinates 0i,02, as illustrated 
in Figure 4. Updating the lagged trigonometric coefficients in accordance with the 
formulas presented in Figures 3a-3e produces trigonometric functions Vi{6) and V2{6) 
which give the values of the ongoing approximations at 9i and ^2, respectively, and 
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whose coefficients are e icit inhomogeneous linear functions of ei eTA> namely, 

Vi{9i) = vie + viu euA + vit Cta, for i = 1, 2. 

The coefficients are explicitly described most concisely in pseudo-code, where e = e^i, 
alp = a, bet = 0, gam = 7, bpd = b + d and ape = a + c: 

if (a > c){ 

axl=alp+(d*d-c*c+2*(a*c-b*d)+a*a-b*b)*2*e; 

ayl=2*(b*b-a*a); 

bxl =bet+ (c*d-a*d-b*c-a*b) *4*e; 

byl=4*a*b; 

cxl=gam+(2*(a*c+b*d)-c*c-d*d+a*a+b*b)*2*e; 

cyl=-2*(a*a+b*b); 

ax2=alp+4*{d*d-c*c)*e; 

ay2=2*(c*c-d*d); 

bx2=bet+8*e*c*d; 

by2=-4*c*d; 

cx2=gam-4*e*(c*c+d*d); 

cy2=2*(c*c+d*d); 

} . 

else{ 

axi=alp+4*e*(a*a-b*b); 
ayl=2*(b*b-a*a); 
bxl=bet-8*e*a*b; 
byl=4*a*b; 

cxl=gam+4*e* (a*a-f b*b) ; 
cyl=-2*(a*a4-b*b); 

ax2=alp+(a*a-b*b+2*(b*d-a*c)-c*c+d*d)*2*e; 
ay2=2*(c*c-d*d); 

bx2=bet+(a*d-a*b+b*c4-c*d)*4*e; 
by2=-4*c*d; 

cx2=gam+(a*a+b*b-2*(a*c+b*d)-c*c-d*d)*2*e; 

cy2=2*(c*c+d*d); 

} 

vlc=cxl+ (((b+bpd)*(b+bpd)-(a+apc)*(a+apc))*axl 

+2*{b+bpd)*{a+apc)*bxl)/((b+bpd)*(b4-bpd)+(a+apc)*(a+apc)); 

vlt=cyl+(((b+bpd)*(b4-bpd)-(a+apc)*(a+apc))*ayl 

+2*byl*(b4-bpd)*(a4-apc))/((b+bpd)*(b+bpd)+(a+apc)*(a+apc)); 

vlu=8./((b+bpd)*(b+bpd)+(a+apc)*(a-fapc)); 

v2c=cx2+(((d+bpd)*(d-hbpd)-(c+apc)*(c+apc))*ax2 

+2*(d+bpd)*(c+apc)*bx2)/((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

v2u=cy2+(({d-l-bpd)*(d-fbpd)-(c+apc)*(c+apc))*ay2 

4-2*(d+bpd)*(c+apc)*by2)/((d+bpd)*(d+bpd)4-(c+apc)*(c+apc)); 

v2t=-8./((d+bpd)*(d+bpd)+(c+apc)*(c+apc)); 

It is easy to minimize [/(^i) - Vi{ei)f + [/(^a) - ^2(^2)]^ as a function of euA, ^ta 
subject to the given constraint, and this procedure provides a suitable regularization. 

12 
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Weighted least-squares .s or other statistical treatments using i it data values at 
further generations allow alternative regularizations, 

Two other simple regularizations which have also been useful are: to require that 
cxjA = 0 if c > a and erA = 0 if a > c, or, to require that cua/c-ta = —[a{a-{-c) + h{h+ 
d)]f[c{a + c)+d[b + d)]; the latter regiilarization is related to imposing differentiability 
of the approximation at Oq. 

The final technical point involves the stopping criteria and terminal processing for the 
cascade of arrow-structures. The basic stopping parameter is NVAN, where a branch 
of the cascade terminates whenever there have been NVAN consecutive generations 
of offspring whose contributions to all coefScients Cn in the bandwidth N have been 
of magnitude at most EPS. There are further parameters governing stopping criteria 
which may be introduced depending upon the scale and resolution of features in the 
input data, including a minimum or maximum number of generations, the size of the 
angle subtended by the chord of a terminal arrow-structure, and so on. 

Another useful stopping criterion is to require that IcaI^H^^aH be negligible for sev- 
eral generations, where Ij^^^ijl denotes the L^-norm of -Qa for which there is the a 
priori estimate W'OaW < lOjac + bd\~^^'^\ that is, one terminates the cascade when 
NVAN consecutive generations of coefficients ca have satisfied the condition that 
|e^|<iEPS \ac+hd\'>/^. 

When NVAN generations of oflEspring contribute negligibly to the specified bandwidth, 
the cascade is terminated with those offspring of greatest generation. Suppose that A 
labels the arrow of such a terminal offspring, adopt the notation in Figure 4, and let 
J C. C denote the circular arc with endpoints , g^ff which contains none of the 
points —1, — i, 4-1 6 C. The final update Tj = alpi cos 6 + beti sin 6 + gami, for i = 
1, 2, of the running trigonometric function on J in accordance with Figures 3a-3e is 
prescribed in pseudo-code with notation as before as follows: 

if(a > c){ 

alpl=alp+e*2*{a*(a+c)-l-d*(d-b)-b*(b+d)-c*(c-a}); 

bet l=bet-he*4* ( c* (d-b)-a* (d-l-b)) ; 

gaml=gam4-e*2*(a*(a-f-c)+c*(a-c)-|-b*(b-f-d)-d*(d-b)); 

alp2=alp4-e*4*(d-c)*(d-f-c); ^ 

bet2=bet-he*8*c*d; 

gam2=gam-e*4*(c*c+d*d); 

} 

else{ 

alpl=alp-fe*4* (a-b) * (a-fb) ; 

betl=bet-e*8*a*b; 

gaml=gam4-e*4*(a*a-f-b*b); 

alp2=aip4-e*2*(a*(a-c)-hb*(d-b)-c*(c+a)+d*(b-^d)); 
bet2=bet+e*4* (a* (d-b)+c* (b+d)) ; 
gam2=gam-l-e*2*(a*(a-c)-c*(a4-c)-b*(d-b)-d*(b-Hd)); 
} 

Now serially letting r denote the trigonometric functions ri and to, again adopt the 
notation of Figure 4, where A labels the arrow corresponding to r. The updated r 
is compared with another trigonometric function cr defined by extrapolation which is 
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uniquely determined b:. amanding that it agree with / at the thr* angles comprising 
the next two generations, B = 9q,Bi,92 in the notation of Figure 4. The Fourier 
coefficients of the piecewise trigonometric function which takes values o{d) — r{Q) on 
the circular arc J and value zero otherwise are easily computed using standard formulas 
and added to those of the ongoing approximation. 

Other algebraic or statistical samplings of higher generation offspring allow alternative 
final processings. 

For instance, one might save computational expense and simply terminate the cascade 
with no final processing at all. 

Furthermore, the terminal arrow-structures could be stored for restart or iterative 
refinement capabilities. 

One might alternatively record in an arrow-structure different transforms of the trigono- 
metric function defined by ct,/3,7. For instance, it is useful to take advantage of the 
renormalization property of wavelets and change these coefiicients by transforming the 
vector field [a cos fi' + Z? sin ^+7]^ by the change of coordinates M^^. This technique 
avoids the necessity of calculating certain large integers attendant to the calculation of 
wavelet coefficients as will be illustrated later. 

One favorable aspect of the method is its advantageous mix of floating-point and integer 
operation types. Moreover, except for the coefficients co, c±i, the implementation of the 
algorithm is purely algebraic, that is, requires only addition and multiplication. Fur- 
thermore, by its very nature as a cascade, the algorithm is amenable to parallelization 
and efficient hardware implementation. 

This completes the description of the implementation of the method of calculating 
Fourier coefficients using arithmetic wavelets. 

Several points will next be addressed which are special to the implementation of the 
methods using modular wavelets. Again there is a basic Step 1, the calculation of 
the modular wavelet transform where the input function is approximated as f{6) « 
1Z/iPa^/i(^)i and a basic Step 2, substitution of known expressions for the Fourier 
coefficients of the modular wavelets; these two steps are again conveniently merged 
into a binary cascade of arrow-structures in practice. 

To elucidate Step 1 for modular wavelets, refer again to Figure 4 for the notation near 
the arithmetic arrow labeled by the matrix A. The single value of the input func- 
tion /(^o) together with the ongoing calculation of the updated trigonometric function 
Or cos 6 + 13 sin 5 -f- 7 in the arrow-structure this time uniquely determines the modular 
wavelet coefficient 

^ S{6o) - {a cos $0+13 sin Bp + 7) 

since V'b(^o) = 0 for every labeling matrix B of generation greater than A; there is 
thus no need for regularization with modular wavelets. 

As to Step 2, recall that a basis of modular wavelets was given before by a collection 
iPb{0] of wavelets where A is the label on an arithmetic arrow and either B = A ov 
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B — SA. In each case, i can calculate that tpsiG) has Fourier e: nsion 



2 (c2 + d2) 



2 (a^ + &2) 




Thus, the method of calculating Fourier transforms using modular wavelets is easily 
derived from the method using arithmetic wavelets. In fact, it is a simple matter to in- 
corporate the two minor modifications just described and alter the included source code 
employing arithmetic wavelets to produce source code employing modular wavelets. 

In the same way for fan wavelets, there is again a Step 1 and Step 2, which are merged 
into a binary cascade of arrow-structures. As for arithmetic wavelets, Step 1 for fan 
wavelets again requires a regularization scheme. As to Step 2, the formulas for Fourier 
coefficients of fan wavelets can be derived from those given before for modular wavelets 
using the identity (^)a{&) = i>A{0) - i^UAid) mentioned before. Again, the included 
source code is easily modified to employ fan wavelets. 



Flow Charts for the Fourier Transform 

Flow charts for various procedures used in the calculation of Fourier transforms are 
presented in Figures 5-8. For definiteness, the flow charts and source code refer to 
arithmetic wavelets without renormalization, but the renormalizing source code will 
also be included for completeness. For purposes of brevity in these flow charts, the 
arrow-structures defined before are referred to in these figures simply as "arrows" . 

In Figure 5 is presented a flow chart for the main driving routine. The input data is 
normalized as indicated in program segment 1. There are separate recursions estab- 
lished in program segments 2 and 3 for the top and bottom of the circle respectively. 



15 



wo 00/44104 



PCT/US99/30584 



Each recursion is estab'. .ed with a call to the subroutine gener, w' h'then calls itself 
in turn. The Fourier coefficients are stored internally with an overall suppressed factor 
of TT, and output data normalization of multiplying by ir is accomplished in program 
segment 4. Output Fourier coefficients are finally displayed before exiting in program 
segment 5. Program segments 2 and 3 are entirely independent and could be performed 
in parallel; more generallj', each of program segments 2 and 3 could be decomposed 
further into multiple parallel procedures. 

A flow chart for the main recursive routine gener is given in Figure 6. The procedure 
starts with a test in program segment 6 to determine if: 

• the argument envy >= NVAN, in which case there have been NVAN 
consecutive generations of negligible contributions, and 

• the argument generation >= MING, in which case there have been a 
required minimum number of iterations in the recursion. 

If both of these inequalities hold, then the procedure passes to program segment 15, 
where the cascade is terminated and the output data corrected with a call to the 
subroutine prune. 

In the contrary case that the recursion should continue, the procedure passes to program 
segment 7, where the descendant arrows in the cascade are determined using the least- 
squares fit to the next generation of data as described before to compute the regularized 
wavelet coefficients of the descendants. These are combined with integer calculations 
to update the lagged trigonometric functions. The procedure then passes to program 
segment 8, where the ongoing approximations to Fourier coefficients are updated to 
include contributions from the two descendant arrows with a call of the subroutine 
Charles for each descendant. The procedure continues with a test in program segment 
9. If the contribution calculated in the subroutine charles for either descendant to any 
Fourier coefficient in the bandwidth was non-negligible, then the recursive argument 
envy is set to zero in program segment 10. In the contrary case that both contributions 
to all Fourier coefficients in the bandwidth were negligible, the control parameter envy 
is increased by one in program segment 13. 

In any case, the procedure passes to program segment 11, where there is a test on the 
number of generations. In case too many iterations of the recursion have occurred, 
it is terminated in program segment 12, and suitable warning of non-convergence of 
the method is written to the output. In the contrary case that the maximum number 
of generations has not been reached, program segment 14 establishes the recursion by 
calhng gener once for each of the descendant arrows with the updated control parameter 
envy and an incremented generation. 

In Figure 7 is presented a flow chart for the subroutine charles. Calling the subroutine 
chaxles has the effect of updating the ongoing approximations to the Fourier coefficients 
for a single argument arrow and returning a flag which keeps track of whether any 
such contribution has been non- negligible. The flag is cleared in program segment 16, 
and several preliminary calculations are accomplished in program segment 17. The 
procedure passes to program segment 18, where it is determined whether to calculate 
the 0,+l,-l Fourier coefficients. If these are to be calculated, then the procedure 
passes to program segment 19, which calls a subroutine to update these three Fourier 
coefficients. Program segments 20 and 21 set the flag as required depending upon 
whether these three contributions are non-negligible. 
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In any case, the proce re passes to program segment 22. If the jp over the band- 
width is complete, then all Fourier coefficients have been updated, so the subroutine 
Charles returns the flag in program segment 23. In the contrary case, the current Fourier 
coefficient is updated in program segment 24. The flag is set as required depending on 
whether the current contribution was non-negligible in program segments 25 and 26. 
The iteration over bandwidth is established by the update in program segment 27 and 
the return to program segment 22. 

In Figure 8 is presented a flow chart for the subroutine prune, which modifies the ar- 
ray of Fourier coefficients when terminating the cascade. The final update of lagged 
trigonometric functions to produce ri,T2 is performed in program segment '28. There is 
one such function for each possible descendant of the argument arrow. The procedure 
passes to program segment 29, where further input data is sampled in order to com- 
pute two trigonometric extrapolations ai,a2, one such extrapolation for each possible 
descendant of the argument arrow. Each function ai — for i — 1,2, is truncated as 
described before, and the Fourier coefficients of the truncated functions are added to 
the ongoing approximations of Foiirier coefficients in program segment 30. Program 
segments 31 and 32 implement the calculation of 0,-}-l,-l Fourier coefficients if desired, 
and the subroutine prune is terminated with the return in program segment 33. 



The Inverse Fourier Transform 

Included in the section entitled "Source Codes" is an implementation in the computer 
language C of a preferred embodiment of the method described in this section employing 
arithmetic wavelets. 

Since the method for calculating the inverse Fourier transform is similar in spirit and 
execution to that for the Fourier transform described in detail before, it need only be 
briefly discussed. The input data includes the specification of a collection of Fourier 
coefficients Cn in a given bandwidth iV. 

There are again two main steps in the calculation: 

Stepl : Calculate arithmetic wavelet coefficients from the given Fourier 
coefficients c^. 

Step 2: Use the wavelet coefficients from Step 1 and the reconstruction al- 
gorithm to output values of the function Y1,a ^a'&a{&) at the Farey quadra- 
ture points on the circle. 

For Step 1, there is again an explicit formula from previous theoretical work. Deep 
mathematical results prove that for n 7^ 0, ±1, 

e'"^ = cos ne -h i sin nd 
where the sum is over all arithmetic arrows, the underlying chord of which has endpoints 
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^, 77 e C, and 



-1, 


n = 0(4) 


0, 


n = 1(4) 


-1, 


n = 2(4) 


0, 


n = 3(4) 



n = 0(4): 
n = 1(4): 
n^2(4): 
n = 3(4) 



Substituting this expression for e' 
gives 



into the given Fourier expansion f(6) 



E 0(4): 
= 1(4): 
= 2(4) 
i 3(4). 



4 E ^f(r+,-)-f-^(r 



for ^ = 0,±1. 



This expression f{6) « 4 + cje^^ + c'.ie-'^ + t?A(^) of f{e) constitutes Step 1 

in the calculation of inverse Fourier transform. 

In fact, the calculation of wavelet coefficients Ba in Step 1 can be implemented using 
mostly integer operations to improve run-time dramatically. To see this, notice that 
the formulas ^ = |±|f , tj = ^ for the endpoints of the chord underlying the arrow 

labeled by the matrix A — may be substituted into the expression for ba to 

yield 



„ _ [{bd - ac) +i{ad-¥ 6c)]" 

(a2+62)n(c2 + ii2)n ^ 

\n{{hd + ac + i)" + ihd + ac- ?)"} + i{bd + ac){{bd + ac + i)'^ - {bd + ac- i)^} j . 

In the obvious implementation of this formula to calculate using integer operations, 
the intermediaxy integer expressions grow too large to be practical for n large if the 
arrow corresponding to A is of large generation. 

This difficulty is overcome by wavelet renormalization, where these attendant large 
integer calculations are obviated entirely as is illustrated in the included source code. 

Step 2 in the method for calculating inverse Fourier transforms depends upon the 
reconstruction algorithm to output function values _on the circle taken by the linear 
combination f{e) ^ c'q + c\e^^ + c!_-^e-^^ + YLa '^a{6). The finiteness property of 
the reconstruction algorithm allows the exact calculation of these function values at 
the sample points Q C C in their ordering determined by the Farey quadrature. The 
two steps are again conveniently merged into a single binary cascade as follows. 
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There is only one con i parameter, which is called SCALE, v. re the method is 
required to determine at least one output value in each circular arc subtending an 
angle SCALE. Thus, as SCALE is decreased the grid of output values on the circle is 
refined. A binary cascade is established using arrow-structures which is similar to that 
for the Fourier transform before. However, the wavelet coefEcients are given in this case 
in closed form by the formula for ca (so no least-squares fit is required). Because of the 
finiteness property of the reconstruction algorithm, exact output values are determined 
at each stage of the recursion. The cascade is terminated with an arrow whenever the 
chords underlying both descendant arrows subtend angles less than SCALE. Output 
data can be written to an array either when calculated to store data in the ordering of 
the Farey quadrature or when terminating the cascade to store data in the standard 
ordering on the circle. Furthermore, the terminal arrow-structures could be stored for 
restart or iterative refinement capabihties. 

This completes the description of the implementation of the method of calculating the 
inverse Fourier transform using arithmetic wavelets. 

Several points will next be addressed which allow the extension of the methods to 
modular and fan wavelets. Again there is a basic Step 1, the calculation of wavelet 
coefficients from Fourier coefficients, and a Step 2, the reconstruction algorithm; in 
each case of .fan or modular wavelets, these two steps are again conveniently merged 
into a binary cascade of arrow-structures in practice. The finiteness condition on fan 
or modular wavelets mentioned before renders Step 2 entirely analogous to that for 
arithmetic wavelets. As to Step 1, simply substitute the identities "Oa = i>A — <f>UA = 
i'UA — ^"^A + described before into the expression above for e'"^ in terms of 

arithmetic wavelets to immediately derive corresponding expressions in terms of fan 
and modular wavelets. 

Thus, the method of calculating inverse Fourier tremsforms using fan or modular 
wavelets is easily derived from the method using arithmetic wavelets. In fact, it is 
a simple matter to incorporate the minor modifications just described and alter the 
included source code employing arithmetic wavelets to produce source code employing 
fan or modular wavelets. 



Flow Charts for the Inverse Fourier Transform 

Flow charts for various procedures used in the calculation of inverse Fourier transforms 
are presented in Figures 9-10. For definiteness, the flow charts and source code refer to 
arithmetic wavelets without renormalization, but the re-normalizing source code will 
also be included for completeness. For purposes of brevity in these flow charts, the 
arrow-structures defined before are referred to in these figures simply as "arrows" . 

In Figure 9 is given a flow chart for the main driving routine. Modified Fourier coef- 
ficients Cq,c[,c'_i are computed in program segment 34. The formula given in Step 1 
for the inverse Fourier transform is applied in program segment 35 to calculate the 
arithmetic wavelet coefficients in terms of Fourier coefficients for a particular family 
of nine arrows. These nine arithmetic wavelet coefficients are required to initialize the 
trigonometric functions for recursions established in program segment 36 with calls 
to the recursive subroutine gener. Recursions with initial conditions corresponding to 
A = U,T,T-^, C/-i,T-\T-^D-\ are undertaken in program segment 36. These 
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recursions, as well as ft .ler possible subrecursions derived from m, are indepen- 
dent and may be undertaken in parallel. This elaborate driving routine is convenient 
because e'"^ is expressed in Step 1 of_the inverse Fourier transform method in terms 
of the normalized arithmetic wavelets which differ from the arithmetic wavelets -Oa 
themselves only for A = C/-^, y-i. 

The procedure passes to program segment 37 to perform an overall correction of the 
function values using the modified Fourier coefficients Cq,c[,c'_i. Output data is dis- 
played and the procedure terminated in program segment 38. 

In Figure 10 is given a flow chart for the main recursive routine gener. In the cas- 
cade, the two descendant arrows of the argument arrow are generated in program 
segment 39 employing Step 1 of the inverse Fourier transform method to calculate the 
wavelet coefficients of the descendants. These expressions are then used to update the 
lagged trigonometric functions in program segment 39. A test is performed in program 
segment 40 to determine if the chord underlying each descendant arrow subtends a 
sufficiently small angle, as estimated by the comparison of an integer quantity with 
SCAT; the relation with the control parameter SCALE discussed before is given by 
SCAT = [exp sinh~^(l/SCALE)]2. In case one of the chords subtends too large an 
angle, the procedure passes to program segment 41. The recursion is established in 
program segment 41 with two calls to gener, where the arguments are given by the 
two current descendant arrows. In the contrary case, the procedure performs the final 
update of the lagged trigonometric functions in program segment 42, then stuffs the 
output array with the new function values in program segment 43, and finally returns 
in program segment 44. 



Data Compression 

As with any method for data compression based on transform coding, five basic ma- 
nipulation are required, namely: 

• wavelet filtering, 

• quantization, 

• storage and retrieval or transmission 

• de-quantization, and 

• reconstruction. 

The wavelet filter has already been fully disclosed as Step 1 for the calculation of the 
Fourier transform, and the wavelet inverse filter or reconstruction algorithm has likewise 
already been described as Step 2 for the calculation of the inverse Fourier transform. 
Application-specific quantization is done according to psychovisual or psychoacoustic 
thresholds. Quantization and storage are furthermore merged with wavelet filtering 
into a single binary cascade as before, and retrieval or transmission are likewise merged 
with reconstruction into another single binary cascade. The reconstruction algorithm 
is exact. Furthermore, source code for it may be transmitted along with compressed 
data since the coding is sufficiently abbreviated and low-level. 

As is well-known, useful compression by transform coding requires an efficient specifi- 
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cation of the basis eier. .ts. This is especially advantageous for new wavelets: a 
top arrow of generation g labeled by the two-by-two integral matrix A is specified by 
g bits, namely, by writing A uniquely as the matrix product of factors U or T. For 

example, the matrix = ^ jq 7^ factors as 

a??)=a 0(3 o(o 00 ?)■ 

so the coefficient for the wavelet is coded by the binary sequence TUUUTTU. 

The method is also especially well-suited to progressive picture build-up or other it- 
erative refinement as follows. The Farey quadrature determines a linear ordering 
on the set Q of quadrature points in the circle C as before. The ordered subset 
Q\{—i, ±1} C Q C C7 is put into bijective correspondence with the collection of all 
arithmetic arrows (where \ denotes the set-theoretic relative complement); as illustrated 

in Figures la-lc, the arithmetic arrow labeled by ^ = corresponds to the com- 

plex number 

g±|±ig±£iec, for a top arrow; 

(bZdj-i(ale) ^ ^' ^ bottom arrow; 
i e C, for the doe. 

Thus, the ordering from the Farey quadrature on Q\{—i, ±1} determines a correspond- 
ing ordering on the set of all arithmetic wavelets. Storage, retrieval, transmission, and 
reconstruction of wavelet coefficients is accomplished at a specified input or output data 
resolution in this canonical ordering since energy compacts to the lesser wavelet coef- 
ficients. In the case of real-time compression and transmission or other manipulation 
not requiring retrieval, the method can be further improved by the manipulation of 
previously computed wavelet coefficients in parallel with the calculation of subsequent 



Complex and Multi-Dimensional Data 

There is nothing remarkable about the extension of the invention to complex-valued 
input data f(9) = u{9) -f- i v(d), where u and v are real- valued functions on the circle 
C. Concentrating for definiteness on arithmetic wavelets, one simple approach is to 
separately and independently transform u and v, 

u{e) ^Y.^A ^A{e) and v{e) « X^/a ^A{e), 

using the techniques already described in order to calculate 

Another direct approach is to allow complex values for the input function / and solve 
for complex wavelet coefiacients using the same algebraic formulas as before. 
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rurning now to inulti-> lensional input/output data, suppose th there are R > 1 
input functions each of which lias M > 1 independent variables, and define the torus 
C^^ to be the iV/-foId Cartesian product 

M times 

Coordinates on C^'^ are given by M-tuples ^= (^i, ^2. • • • , Om) of angles, and the input 
data is specified by an i?-tuple of multivariable real- or complex- valued functions 

F{e)=={f\o)j''{e):---J''m 

Using the bijection described in the previous section on data compression, each quadra- 
ture point Z G <5\{— i, ±1} uniquely determines a corresponding arithmetic wavelet 
'dzi^y, by convention, one also defines = i5±i(5) = 1. Concentrating on arith- 

metic wavelets for definiteness, define the corresponding multiwavelet 

As before, the input function F must be normalized. To this end, define a trigonometric 
monomial to be a product of the form X'^{6i)X'^{62)--- X^{9m), where each X*(0) is 
given by one of X^{e) = cos 6, X^{$) - sin 0, or X^{0) = 1. The values of F at 
the points in {— i,±l}^ uniquely determine coefficients in a linear combination u(0) 
of trigonometric monomials so that 1/(6) — F{$) for each point of {-i,±l}^. Just 
as in the one-dimensional wavelet filter already discussed, now one must consider the 
normalization F(e) = F{e) - u{e), so F is zero at each point of {—i, ±1}^, 

Adopting the Farey quadrature in each factor circle C of the torus , there is an 
inducedjexicographic ordering on the set Q^^\{-z, ±1}^ C of potential sample 
points Z € C^^ . One separately approximates each coordinate function 

f{e) « X! 4 ^^(i^)' for 2 = 1,2,...,^, 

of F{9) as a finite linear combination of multiwavelets with application-specific zonal 
sampling for optimum energy compaction; this calculation may employ least-squares or 
other regularization in each factor as well as renormalization as for the one-dimensional 
wavelet filter. The methods are again elegantly implemented in practice as M nested 
binary cascades. The final processing for the cascades of arrow-structures in the cal- 
culation of multi-dimensional Fourier transforms may involve adding suitable linear 
combinations of trigonometric monomials based on further sampling, or one may sim- 
ply truncate with no final processing at all, as in the one- dimensional case. 

These same remarks apply verbatim to fan and modular wavelets. 



Mathematical Basis 

"The decorated Teichmiiller space of punctured surfaces" , Communications in Mathe- 
matical Physics 13, pp. 299-339 (1987), written by the author of this patent applica- 
tion, began the study of certain abstract geometric coordinates called lambda lengths 
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in the context of hype Aic geometry on Riemann surfaces. Th» leas developed in 
this and other related papers have found applications in high energy physics. A more 
recent publication in this series of about 20 papers is entitled "Universal constructions 
in Teichmiiller Tlieory", Advances in Mathematics 98, pp. 143-215 (1993). Here it 
was shown how to generalize lambda lengths to provide coordinates on a certain space 
of homeomorphisms (i.e. continuous bijections whose inverses are also continuous) of 
the circle, where there is one generalized lambda length for each arithmetic arrow. 
In further recent work described in "The Lie algebra of homeomorphisms of the cir- 
cle", Advances in Mathematics 140, pp. 282-322 (1998), written jointly with Feodor 
Malikov, there is described a representation of the coordinate deformations of these 
generalized lambda lengths as explicit vector fields on the circle. These vector fields 
are called "normalized elementary vector fields", and their study led to other vector 
fields on the circle called "fans" and "hyperfans" . These three families of vector fields 
correspond, respectively, to the arithmetic, fan, and modular wavelets described here 
via the identification of the vector field f{0)-§Q with the function f{d); there was, how- 
ever, no discussion of wavelets, sampling of data, approximation of functions, or any 
algorithms in the published literature. 

Each arithmetic wavelet is once-continuously differentiable on the circle, compactly 
supported, and locaUzed in space. Arithmetic wavelets are not compactly supported 
in frequency, but instead, the frequency profile of an arithmetic wavelet is given alge- 
braically by the formula for used in Step 2 in the calculation of Fourier transforms; 
a non-compactly supported locaHzation in frequency follows directly from this. The 
asserted formula for can be derived without much difficulty but in several cases 
directly from Figures 3a-3e integrating twice by parts the usual expression for Fourier 
coefiicients using the fact is once-continuously differentiable. More difiicult to derive 
is the formula given for the exponential functions e'"^ in terms of arithmetic wavelets 
which was used in Step 1 of the calculation of inverse Fourier transforms; this subtle 
computation with lambda lengths is given in "The Lie algebra of homeomorphisms of 
the circle". 

Similar remarks apply to fan and modular wavelets, but fan wavelets are only contin- 
uous (they are not differentiable), and modular wavelets are not even continuous. 

In order to explain the failure of orthonormaiity of arithmetic wavelets and illuminate 
the regularization used in the arithmetic wavelet filter, fix some matrix A labeling 
an arithmetic arrow whose underlying chord has initial point q in the circle C. The 
sequence V^A of matrices label the arithmetic arrows whose underlying chords also 
have q as initial point, where n denotes an arbitrary integer. There is one infinite 
linear dependence 

0 = ^^yn^. 



for each q e Q, and this explains the additive degree of freedom in the calculation of 
euA and cta in the wavelet filter which is handled by regularization as was discussed 
before. This additive degree of freedom is also manifest in the initial specification 
6/ = 0 of wavelet coefficient on the doe in the wavelet filter. 

In order to better understand the Farey quadrature, it is convenient to transform the 
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disk D to the upper h plane U = {u + iv: v > Q] hy means of'' transform 

2-1" 

This function i\ is a homeomorphism on the interior of D which maps the unit circle 
C bijectively to the real-axis + : f = 0} together with an additional point oo 
at infinity, where K{-1) = 0, K{-i) = 1, and = oo. Furthermore, as q e Q 

varies over all quadrature points, the real number K{q) varies over all rational numbers, 
where the rational number 2 corresponds to the quadrature point 2^ e Q. In order 
to illustrate the Farey tesselation itself in U, whenever 91,92 S Q are the endpoints 
of the chord underlying an arithmetic arrow in D, draw the semi-circle in U with real 
endpoints 51,92; in the degenerate case that one of the quadrature points is infinite, 
say the point 92 = 00, then draw the vertical ray in U with endpoint Qi . Applying this 
procedure to each arithmetic arrow produces the classical model of the Farey tesselation 
in U that is indicated in Figure 11. This figure illustrates the relationship between 
the Farey tesselation and the indicated circle packing, where the circle in the figure 
that is tangent to the real-axis at the rational point p/g has diameter g~^. For further 
information on number-theoretic aspects of this classical figure, see, for instance. Ford's 
book "Automorphic Functions" , Chelsea (1972). 



Source Codes 

In order to adequately disclose the methods, it is the purpose of this paragraph to 
include source codes for their implementations in the computer language C. Source 
code is presented for each of the following two implementations. 

• awft is the implementation of a preferred embodiment of the method 
for computing the Fourier transform of real-valued input data defined on 
the circle. The algorithm depends upon a bandwidth N, tolerance EPS, 
generation cut-off NVAN, and minimum generation MING as described 
before. A further flag SMODE=l determines that the output data will 
include the 0,4-1,-1 Fourier coefficients, and SMODE=0 otherwise. (The 
only configuration in either code which requires any library routines is 
SM0DE=1, which requires math.h.) 

• awift is the implementation of a preferred embodiment of the method 
for computing the inverse Fourier transform of specified complex- valued 
Fourier coefficients defined in a specified bandwidth N. The output data 
is controlled by a single parameter SCAT, which is related to the param- 
eter SCALE discussed before by SCAT=[exp sinh * (-1)(1/SCALE)] ' 2. 
(SCAT varies inversely with SCALE and may be compared with conve- 
niently derived approximate quantities.) 

The comparison of the subroutines tgener in awft and in awift illustrates the easy 
transition between source codes for real and for complex data. The corresponding 
coraplex-awft and real-awift are therefore easily derived from the included source codes. 
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Furthermore, the met 1 for compression of one-dimensional in] data amounts to 
the first step of awft, then quantization/de-quantization, then the second step of awift. 
The included source codes together thus also implement the method for compression 
of one-dimensional data since source code for it may be extracted from the included 
source codes. Finally, it is straight-forward to write the further driving routine required 
for multi-dimensional data compression, so the source codes included with this patent 
application are sufficient to readily implement the method for multi-dimensional data 
compression as well. 

The included implementations have not been optimized. (For instance, straight-forward 
optimization has improved run-time by as much as 70% over the included awft code; 
on the other hand, optimization destroys intelligibility of the coding in relation to 
the underlying algorithm, and this explains the inclusion of unoptimized source code.) 
Another substantial improvement in run-time over the included source codes can be 
achieved by factoring so as to replace various manipulations of complex numbers with 
manipulations of integers as discussed before. 

Moreover, for clarity, the included source codes do not take advantage of renormal- 
ization, but for completeness, to each of awft and awift is appended source code to 
replace the subroutines tgener and bgener by corresponding subroutines which do em- 
ploy renormalization. 
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